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An analytical description of the interface motion of a collapsing nanometer-sized spherical cavity 
in water is presented by a modification of the Rayleigh-Plesset equation in conjunction with ex- 
plicit solvent molecular dynamics simulations. Quantitative agreement is found between the two 
approaches for the time-dependent cavity radius R(t) at different solvent conditions while in the con- 
tinuum picture the solvent viscosity has to be corrected for curvature effects. The typical magnitude 
of the interface or collapse velocity is found to be given by the ratio of surface tension and fluid vis- 
cosity, v ~ 7/77, while the curvature correction accelerates collapse dynamics on length scales below 
the equilibrium crossover scales (~lnm). The study offers a starting point for an efficient implicit 
modeling of water dynamics in aqueous nanoassembly and protein systems in nonequilibrium. 



I. INTRODUCTION 

Hydrophobic hydration in equilibrium is a phe- 
nomenon which exhibits qualitatively different behavior 
at small and large length scales. 1,2 While small solutes 
(radii R <lnm) are accommodated by water with only 
minor perturbations, larger solutes (i? >lnm) induce ma- 
jor rearrangements of water interfacial structure. As a 
consequence the solvation free energy G(R) of small hy- 
drophobic cavities scales with solute volume while for 
larger cavities it grows with surface area (as a good 
approximation near liquid-vapor coexistence) accompa- 
nied by weak solvent dewetting at extended restrain- 
ing hydrophobic surfaces. 3 Furthermore, water, which 
is close to the liquid-vapor transition at normal condi- 
tions, can minimize interface area by locally evaporat- 
ing and forming a 'nanobubble' within hydrophobic con- 
finement. Evidence of bubble formation in confined ge- 
ometry has been given early by computer simulations 
of smooth plate-like solutes, 4 but more recently it has 
been demonstrated in varying degrees in atomistically 
resolved plate-like solutes,— £ hydrophobic tubes and ion 
channels, 7 ' - and in the collapse of proteinsj^jlfi, suggesting 
that it plays a key role in the stabilization and folding 
dynamics of certain classes of biomolecules j 1 1 1 12 Experi- 
mental evidence of nanobubbles in strong confinement (in 
contrast to bubbles at a single planar surface 3 -) has been 
given for instance in studies of water between hydropho- 
bic surfaces ji 3 - in zeolites and silica nanotubeS ) 14 i 15 and 
on a subnanometer scale in nonpolar protein cavities^ 

The dewetting induced change in solvation energy is 
typically estimated using simple macroscopic arguments 
as known from capillarity theory, e.g. by describing in- 
terfaces with Laplace- Young (LY) type of equationsi 14 i 17 
Recently an extension of the LY equation has become 
available which extrapolates to microscopic scales by in- 
cluding a curvature correction to the interface tension 
and considering atomistic dispersion and electrostatic 
potentials of the solvated solute explicitly^ Although 
those macroscopic considerations (e.g,. the concept of 
surface tension) are supposed to break down on atom- 
istic scales they show surprisingly good results for the 



solvation energy of microscopic solutes, e.g. alkanes and 
noble gases, and quantitatively account for dewetting 
effects in nanometer-sized hydrophobic confinement. 19 
While we conclude that the equilibrium location of the 
solute-solvent interface seems to be well described by 
those techniques, nothing is known about the interface 
dynamics of evolution and relaxation. In this study we 
address two fundamental questions: First, what are the 
equations which govern the interface motion on atomistic 
(~lnm) scales? Secondly, does the dynamics exhibit any 
signatures of the length scale crossover found in equilib- 
rium? 

On macroscopic scales the collapse dynamics of a (va- 
por or gas) bubble is related to the well-known phe- 
nomenon of sonoluminescence^ The governing equa- 
tions can be derived from Navier-Stokes and capillarity 
theory and are expressed by the Rayleigh-Plesset (RP) 
equation^ We will show that the RP equation simpli- 
fies in the limit of microscopic cavities and can be ex- 
tended to give a quantitative description of cavity inter- 
face dynamics on nanometer length scales. We find a 
qualitatively different dynamics than the typical "mean- 
curvature flow" description of moving interfaces, 22 in par- 
ticular a typical magnitude of interface or collapse veloc- 
ity given by the ratio of surface tension and fluid viscosity, 
v ~ 7/77- Our study is restricted to the generic case of 
the collapse of a spherical cavity and is complemented 
by explicit solvent molecular dynamics (MD) computer 
simulations. We note here that recently, Lugli and Zer- 
betto studied nanobubble collapse in ionic solutions by 
MD simulations on similar length scales.- 2 ^ While their 
MD data compares favorably with our results their in- 
terpretation and conclusions in terms of the RP equa- 
tion are different. We will resume this discussion in the 
conclusion section. 

In this study we show that a simple analytical approach 
quantitatively describes microscopic cavity collapse for 
a variety of different solvent situations while the sim- 
ulations suggest that the solvent viscosity needs to be 
corrected for curvature effects. Our study might offer a 
simple starting point for an efficient implicit modeling 
of water dynamics in aqueous nanoassembly and protein 
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systems in nonequilibrium. 



II. THEORY 

The Rayleigh-Plesset equation for the time evolution 
of a macroscopic vapor bubble with radius R(t) can be 
written &s£^ 



,>„, ! RR + Ik 2 



= AP + 4^ + ^, (1) 



where p m is the solvent mass density, AP = P — P v the 
difference in liquid and vapor pressures, 77 the dynamic 
viscosity, and 7 the liquid-vapor interface tension. While 
for macroscopic bubble radii the inertial terms (left hand 
side) control the dynamics, for decreasing radii the fric- 
tional and pressure terms (right hand side) grow in rel- 
ative magnitude and eventually dominate, so that com- 
pletely overdamped dynamics can be assumed on atom- 
istic scales: 
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A rough estimate for the threshold radius Rt below which 
friction dominates is given when the Reynolds number 
1Z = vRp m /r\ becomes unity and viscous and inertial 
forces are balanced. With a typical initial interface ve- 
locity of the order of v ~ 7/77 [from R(0) = in cq. fl}] 
we obtain 



Rt = V 2 /(Pml) 



(3) 



which is ~ lOnm for water at normal conditions. Note 
that this threshold value can deviate considerably for a 
fluid different than water and that the viscosity typically 
has a strong temperature (T) dependence which implies 
that Rt can change significantly with T. 

In equilibrium (R = 0) the remaining expression in 
eq. © is the (spherical) LY equation AP + 2-f/R= 0. 
Thus eq. describes a linear relationship between cap- 
illary pressure and interface velocity where R/(Arj) plays 
the role of an interface mobility (inverse friction) ^ In- 
terestingly, the mobility is linear in bubble radius which 
leads to a constant velocity driven by surface tension in- 
dependent of radius (assuming P ~ 0) ; this has to be con- 
trasted to the typically used capillary dynamics which is 
proportional to the local mean curvature oc 1/R. 22 

Generalizations of the LY equation to small scales are 
available by adding a Gaussian curvature term (~ 1/R 2 ) 
as shown by Boruvka and Neumann^; that has been 
demonstrated to be equivalent to a first order curvature 
correction in surface tension, i.e. j(R) = 7oo(l~ <5t/P), 18 
where St is the Tolman length^ and 700 the liquid-vapor 
surface tension for a planar interface (R — 00). The 
Tolman length has a magnitude which is usually of the 
order of the size of a solvent molecule. Furthermore, it 
has been observed experimentally that the viscosity of 



strongly confined water can depend on the particular na- 
ture of the surface/interfaced We conclude that in gen- 
eral one has to anticipate that - analogous to the surface 
tension - the effective interface viscosity obeys a curva- 
ture correction in the limit of small cavities due to water 
restructuring in the first solvent layers at the hydropho- 
bic interface. In the following we make the simple first 
order assumption that the correction enters eq. @ also 
linear in curvature (~ 1/R) yielding 
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where the constant S v { s is the coefficient for the first order 
curvature correction in viscosity and 7700 the macroscopic 
bulk viscosity. Additionally, we define S = 5 V i S — St and 
second order terms in curvature are neglected. We note 
that the choice of the 1 / P-scaling of the viscosity curva- 
ture correction has no direct physical justification and is 
arbitrary. We think however, that a curvature correction 
based on an expansion in orders of mean curvature is the 
simplest and most natural way for such a choice. 

In water at normal conditions the pressure terms in 
(HI are negligible so that for large radii (R 3> S) the in- 
terface velocity is constant and R(t) = i?o — 7oo/ (2?7oo)i. 
This leads to a collapse velocity of about v ~0.4A/ps 
(40m/s) which is 6% of the thermal velocity of water 
v t h = \JikBT /m showing that dissipative heating of the 
system is relatively weak on these scales. A rough esti- 
mate for the dissipation rate can be made by the released 
interfacial energy dG(R,t)/dt ~ d(4TrR(t) 2 'y 00 )/dt = 
—4:Tr'y 2 c R(t)/r] O0 yielding for instance dG(R,t = 0)/dt ~ 
— 35fcsT/ps for a bubble with Rq =2nm. At small 
radii (R ~ S) the solution of ([J| goes as R(t) ~ 
±y/ const — (£700/ 7/00)/; decreasing or increasing the ve- 
locity depending on the sign of S = S vis — St, i.e. the 
acceleration depends on the particular sign and mag- 
nitude of the curvature corrections to surface tension 
and viscosity. For large pressures and radii the first 
term dominates which gives rise to an exponential de- 
cay R(t) ~ exp[— AP/(4i7oo)t]. While extending to small 
scales we have assumed that the time scale of internal in- 
terface dynamics, i.e. hydrogen bond rearrangements,-^ 
is much faster than the one of bubble collapse. 



III. MD SIMULATION 

In order to quantify our analytical predictions we 
complement the theory by MD simulations using ex- 
plicit SPC/E water.— The liquid- vapor surface tension 
of SPC/E water has been measured and agrees with the 
experimental value for a wide range of temperatures." 
For T = 300K and P = lbar we have 7^0 = 72mN/m. 
The Tolman length has been estimated to be St — 
0.9 A from equilibrium measurements of the solvation 
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energy of spherical cavities^ At the same conditions the 
dynamic viscosity of SPC/E water has been found to 
be ryoo = 6.42 • 10~ 4 Pa-s,— —24% smaller than for real 
water. In experiments in nanometer hydrophobic con- 
finement and at interfaces however, the viscosity shows 
deviations from the bulk value but remains comparable. 26 
We proceed by treating the viscosity "q^ as an adjustable 
parameter together with its curvature correction coeffi- 
cient S v is- 

The MD simulations are carried out with the 
DLPOLY2 package^ 2 - using an integration time step of 
2fs. The simulation box is cubic and periodic in all three 
dimensions with a length of L — (61.1 ± 0.2)A in equi- 
librium involving N — 6426 solvent molecules. Electro- 
static interactions are calculated by the smooth-particle 
mesh Ewald summation method. Lennard- Jones inter- 
actions are cut-off and shifted at 9A. Our investigated 
systems are at first equilibrated in the NPT ensemble 
with application of an external spherical potential of the 
form (3V(r) — [A/(r — R' )] 12 and all molecules removed 
with r < R' Q since vapor can safely be neglected on these 
scales. This stabilizes a well-defined spherical bubble of 
radius Rq ~ R' Q + lA. We define the cavity radius by 
the radial location where the water density p(r) drops to 
half of the bulk density po/2. Thirty independent config- 
urations in 20ps intervals are stored and serve as initial 
configurations for the noncquilibrium runs. We employ a 
Nose-Hoover barostat and thermostat with a 0.2ps relax- 
ation time to maintain the solvent at a pressure P and 
a temperature T. Other choices of relaxation times in 
the reasonable range between 0.1 and 0.5ps do not alter 
our results. In the nonequilibrium simulations the con- 
straining potential is switched off and the relaxation to 
equilibrium is averaged over the thirty runs. 

IV. RESULTS 



system 


P/bar 


T/K 


CN a ci/M Q/e 


W(10~ 4 Pa-s) 


I 


1 


300 








5.14 


II 


1 


300 


1.5 





5.94 


III 


1 


277 








8.48 


IV 


2000 


300 








4.56 


V 


1000 


300 








4.72 


VI 


1 


300 





+2 


5.14 



TABLE I: Investigated system parameters: pressure P, tem- 
perature T, and salt (NaCl) concentration c. In system VI 
a fixed ion with charge Q = +2e is placed at the center of 
the collapsing bubble. The viscosity rjoo is a fit-parameter in 
systems I-V (see text). 

We perform simulations of six different systems I- VI 
whose features are summarized in Tab. I and differ in 
thermodynamic parameters T and P (I, III, IV, and V) 
but also inclusion of dispersed salt (II), and the influ- 
ence of a charged particle in the bubble center (VI) are 



considered. Note that the exact value of the crossover 
length scale (however defined) can depend on the detailed 
thermodynamic or solvent condition but remains close to 
lnm. 2 




r/A 



FIG. 1: Interface density profiles p(r)/po for system I are 
plotted vs. the radial distance r from the bubble center for 
different times i/ps=l,5,10,14,17,19,23. Symbols denote MD 
simulation data and lines are fits using 2p(r)/po = erf{[r — 
R(t)]/d} + l. The bubble radius R(t) is defined by the distance 
at which the density is po/2 (dotted line). The inset shows 
the "10-90" thickness r = 1.8124 d of the interface vs. P for 
initial radii Po = 19.83A (pluses) and Po = 25. 6A (crosses). 

System I is at normal conditions (T=300K, P=lbar) 
and consists of pure SPC/E water. Fig. 1 shows the 
observed interface profiles in the nonequilibrium situa- 
tion at different times i/ps=l, 5, 10, 14, 17, 19, and 23 
starting from an initial radius Po = 19.83A. The liquid- 
vapor interface stays relatively sharp in the process of 
relaxation but broadens noticeably for smaller radii. At 
t ~ 23ps the system is completely relaxed to a homoge- 
neous density distribution. The same time scale of bubble 
collapse has been found in explicit water computer simu- 
lations of dewetting in nanometer-sized paraffin plates, 17 
polymersjii and atomistically resolved proteins^i^ 

We find that the interface profiles can be fitted very 
well with a functional form 2p(r)/po = eri{[r — R(t)]/d}+ 
1, where d is a measure of the interface thickness. The 
interface fits are also shown in Fig. 1 together with the 
MD data. The experimentally accessible "10-90" thick- 
ness r of an interface is the thickness over which the 
density changes from O.lpo to 0.9po and is related to the 
parameter d via r = 1.8124 d. While experimental values 
of t for the planar water liquid-vapor interface vary be- 
tween ~ 4 and 8A the measured values for SPC/E water 
in the finite simulation systems are r m =3 to 4A^ We 
find a strongly radius-dependent function t(R) plotted 
in the inset to Fig. 1 for initial radii Po = 19.83A and 
Po = 25. 6A. For P ~ Po the thickness increases during 
the following 5ps from the equilibrium value r ~ 3A to 
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FIG. 2: Time evolution of the cavity radius R(t) for parame- 
ters as defined in systems I- VI. The solution of the modified 
RP equation Q (lines) is plotted vs. MD data (symbols). 
The inset shows the solution of the modified RP equation in- 
cluding inertia terms, cf. lhs of l[T]l. (dashed lines) compared 
to eq. (O for system I with initial radii Ro = 19. 83 A and 
Ro = 10.0A. 

about r ~ 4.5 — 5 A independent of R . While the exact 
equilibrium thickness at t — depends on the particu- 
lar choice of the confining potential V(r) (e.g., a softer 
potential might lead to a broader initial interface) this 
suggest that 4.5-5A is the typical interface thickness 
for a bubble of lnm size. Regarding the slope of the 
curve one might speculate that t(R — > oo) saturates to 
the thickness Too of the measured planar interface for 
Rq — > oo. For R < 10 A the thickness increases twofold 
during the relaxation to equilibrium. This broadening 
might be attributed to increased density fluctuations and 
the structural change of interfacial water in the system 
when crossing from large to small length scales which has 
been shown to happen in equilibrium at ~ lnmi^ 

In Fig. 2 we plot the time evolution of the bubble ra- 
dius R(t) for all investigated systems. Let us first focus 
on the simulation data of system I (circles). As antic- 
ipated the bubble radius decreases initially in a linear 
fashion while for smaller radii (R(t) < 10A) the velocity 
steadily increases. From the best fit of eq. (0| we find 
a viscosity rjoo = 5.14 • 10~ 4 Pa-s and its curvature cor- 
rection coefficient <5 V ; S = 4.4A. Although investigating a 
confined system with large interfaces the viscosity value 
differs only 20% from the SPC/E bulk value. Further- 
more, from our macroscopic point of view the MD data 
show that high curvature decreases the viscosity and the 
latter has to be curvature-corrected with a (positive) co- 
efficient larger than the Tolman length St- If the surface 
tension decreased in a stronger fashion with curvature 
than viscosity the collapse velocity would drop in qual- 
itative disagreement with the simulation. The overall 
behavior of R(t) and the collapse velocity of about ~ 



lA/ps agrees very well with the recent MD data of Lugli 
and Zerbetto, who simulated the collapse of a lnm sized 
bubble in SPC water ,22 

The inset to Fig. 2 shows the solution of eq. J4]) in- 
cluding incrtial terms [left hand side of |T])] to check the 
assumption of overdamped dynamics. While incrtial ef- 
fects are indeed small but not completely negligible for 
an initial radius Rq — 19.83A they basically vanish for 
Rq = lOA. Interestingly, the inertial effects are not visible 
in the MD simulation data at all. We attribute this obser- 
vation to the finite and periodic simulation box which is 
known to suppress long-ranged inertial (hydrodynamic) 
effects^ 

In the following we assume <5 v i s to be independent of 
the other parameters and treat only r/oo as adjustable 
variable. In system II we add 175 salt pairs of sodium 
chloride (NaCl) into the aqueous solution resulting in a 
concentration of c ~1.5M. The ion-SPC/E interaction 
parameters are those used by Bhatt et al£^ who mea- 
sured a linear increase of surface tension with NaCl con- 
centration in agreement with experimental data. While 
this increment for c = 1.5M is about small 2-3%, the vis- 
cosity has been measured experimentally to increase by 
approximately 18% at 298.15Ki2£ Indeed by comparing 
the simulation data to the theory we find a 16% larger 
viscosity 7700 = 5.94 • 10~ 4 Pa-s. A slower collapse velocity 
has been found also in the MD simulations of Lugli and 
Zerbetto in concentrated LiCl and CsCL solutions when 
compared to pure water ^ 

In system III we investigate the effect of lowering the 
temperature by simulating at T = 277K. While only a 5% 
increase of the water surface tension (SPC/E and real 
water) is estimated from available data£2 the viscosity 
depends strongly on temperature: the relative increase 
has been reported to be between 55 — 75% for SPC/E 
water (85% for real water) Inspecting the MD data 
and considering the surface tension increase we find in- 
deed a large decrease in viscosity of 65% with a best-fit 
T)oc = 8.48 • 10 _4 Pa-s. Both systems, II and III, show 
that solvent viscosity has a substantial influence on bub- 
ble dynamics as quantitatively described by our simple 
analytical approach. In systems IV and V we return to 
T = 300K but increase the pressure P by a factor of 
2000 and 1000, respectively. Best fits provide viscosities 
which are around 10% smaller than at normal conditions 
in agreement with the very weak pressure dependence of 
the viscosity found in experiment s 37 ' 38 at T=300K. The 
major contribution to the faster dynamics comes explic- 
itly from the pressure terms in eq. Q. Although mov- 
ing away from liquid- vapor coexistence by increasing the 
pressure up to 2000bar we assume (and verify hereby) 
that the bubble interface tension can still be described 
by 7oo- 

In system VI we investigate the influence of a hy- 
drophilic solute on the bubble interface motion in order 
to make connection to cavitation close to molecular (pro- 
tein) surfaces. As a simple measure we fix a divalent ion 
at the center of the bubble so that we retain spherical 
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symmetry. The ion is modeled by a Lennard-Jones (LJ) 
potential U LJ (r) = 4e[(cr/r) 12 - [<r/rf] with Q = +2e 
point charges and uses the LJ parameters of the SPC/E 
oxygen-oxygen interaction. As demonstrated recently 
the LY equation can be modified to include dispersion 
and electrostatic solute-solvent interactions explicitly, 18 
which extends Q to 
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The last term in ([5]) is the Born electrostatic energy den- 
sity of a central charge Q in a spherical cavity with ra- 
dius R with low dielectric vapor e v = 1 surrounded by 
a high dielectric liquid (1/ej ~ 0). The electric field 
around the ionic charge and the dispersion attracts the 
surrounding dipolar water what accelerates and eventu- 
ally completely governs the bubble collapse below a ra- 
dius R(t) < 13A (t > 7ps) as also shown in Fig. 2. 
The theoretical prediction (|5|) agrees very well without 
any fitting using the viscosity from system I. We find 
that the acceleration is mainly due to the electrostatic 
attraction; the dispersion term plays just a minor role 
while the excluded volume repulsion eventually deter- 
mines the final (equilibrium) radius of the interface with 
R(t = 00) ~ 2A. 



V. CONCLUSIONS 

In conclusion, we have presented a simple analytical 
and quantitative description of the interface motion of 
a microscopic cavity by modifying the macroscopic RP 
equation. Based on our MD data we find for the macro- 
scopic description that analogous to the surface tension 
the viscosity has to be corrected for curvature effects, a 
prediction compelling to investigate further in detail and 
probably related to the restructuring of interfacial water 
for high curvatures (small R). The viscosity correction 
accelerates collapse dynamics markedly below the equi- 
librium crossover scale (~lnm) in contrast to the pure 
equilibrium picture where surface tension decreases what 



slows down the collapse. Further, we find that the dy- 
namics is curvature-driven due to the corrections to sur- 
face tension and viscosity, not due surface tension as often 
postulated* 2 - 2 - As a simple estimate, the interface velocity 
is typically given by the ratio of surface tension and fluid 
viscosity, v ~ j/rj. 

A comment has to be made regarding the recent work 
of Lugli and Zerbetto on MD simulations of nanobubblc 
collapse in ionic solutions. While their MD data of the 
collapse velocity for a lnm bubble agree very well with 
our results their interpretation in terms of the RP equa- 
tion is different. They fit the 'violent regime' solution 
of the RP equation to the data [which is the solution of 
only the inertial part, left hand side of ((TJ] and argue 
that the violent regime still holds on the nm scale. As 
demonstrated in this work, we arrive to a different con- 
clusion: the collapse is friction dominated, the collapse 
driving force is mainly capillary pressure, and we suggest 
that the microscopic viscosity has to be curvature cor- 
rected to explain the high curvature collapse behavior in 
the MD simulations. The good agreement between our 
modified RP equation and the MD data for different sol- 
vent conditions, leading for instance to an altered solvent 
surface tension or viscosity, support our view. 

We finally note that extensions of the LY equation are 
based on minimizing an appropriate free energy G(R) or 
free energy functiona l 18 ' 24 so that we can write in a more 
general form R ~ [dG(R)/dR]/[r](R)R}. It is highly de- 
sirable to generalize this simple dynamics further to ar- 
bitrary geometries with which a wide field of potential 
applications might open up, i.e. an efficient implicit mod- 
eling of the water interface dynamics in the nonequilib- 
rium process of hydrophobic nanoassembly, protein dock- 
ing and folding, and nanofluidics. 
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